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ABSTRACT 

We present an analytical model for the non-spherical collapse of overdense regions out 
of a Gaussian random field of initial cosmological perturbations. The collapsing region is 
treated as an ellipsoid of constant density, acted upon by the quadrupole tidal shear from 
the surrounding matter. The dynamics of the ellipsoid is set by the ellipsoid self-gravity 
and the external quadrupole shear. Both forces are linear in the coordinates and therefore 
maintain homogeneity of the ellipsoid at all times. The amplitude of the external shear is 
evolved into the non-linear regime in thin spherical shells that are allowed to move only 
radially according to the mass interior to them. The full dynamical equations then reduce to 
a set of nine second-order ordinary differential equations, which reproduce the linear regime 
behaviour but can be evolved past turnaround, well into the non-linear regime. We describe 
how the initial conditions can be drawn in the appropriate correlated way from a random field 
of initial density perturbations. The model is applied to a restricted set of initial conditions 
that are more suitable to the above approximations; most notably we focus on the properties 
of rare high density peaks 2a). By considering many random realizations of the initial 
conditions, we calculate the distribution of shapes and angular momenta acquired by objects 
through the coupling of their quadrupole moment to the tidal shear. The average value of the 
spin parameter, (A) ~ 0.04, is found to be only weakly dependent on the system mass, the 
mean cosmological density, or the initial power spectrum of perturbations, in agreement with 
N-body simulations. For the cold dark matter power spectrum, most objects evolve from 
a quasi-spherical initial state to a pancake or filament and then to complete virialization. 
Low-spin objects tend to be more spherical. The evolution history of shapes is primarily 
induced by the external shear and not by the initial triaxiality of the objects. The statistical 
distribution of the triaxial shapes of collapsing regions can be used to test cosmological 
models against galaxy surveys on large scales. 

Subject headings: cosmology: theory 
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1. INTRODUCTION 



In the standard cosmological model, density perturbations grow by gravitational insta- 
bility from their small initial amplitude to the observed non-linear structure, thus accounting 
for the galaxies, clusters, superclusters, filaments and voids in the present universe (Peebles 
1993). Despite the simplicity of these initial conditions, the analytical understanding of the 
advanced stages of gravitational instability, where the density contrast Sp/p exceeds unity, 
is limited. Numerical codes are frequently used in an attempt to uncover the non-linear 
dynamics of collapsing regions and in order to predict the observational consequences of spe- 
cific models for structure formation. However, simulations are often expensive in computer 
time and are therefore limited in resolution or total volume. An approximate extension of 
linear perturbation theory into the non-linear regime is achieved by the Zel'dovich (1970) 
approximation (see review by Shandarin & Zel'dovich 1989). Another approach which ex- 
trapolates the physics even further, up to the virialization phase of bound objects, is the 
spherical collapse model (Gunn & Gott, 1972; Peebles 1980, §19). Under the assumption of 
sphericity, the non-linear dynamics of a collapsing shell is determined by the mass interior 
to it and described by the parametric solution to the dynamics of a closed universe. Al- 
though the Zel'dovich approximation generically yields pancakes and numerical simulations 
show that many collapsing regions are filamentary (e.g. Park 1990; Bertschinger & Gelb 
1991; Cen & Ostriker 1993), the spherical model became popular because of its simplicity. 
Most notably, it was integrated into the Press-Schechter formalism (Press & Schechter 1974) 
for calculating the mass function of collapsed objects in the universe. In reality, overdense 
regions in galaxy surveys of the local universe (e.g., Geller & Huchra 1989; Maddox et al. 
1990; Saunders et al. 1991; Shectman et al. 1992; Strauss et al. 1992) have complicated 
triaxial shapes with sheets and filaments being common features. Moreover, recent work 
has demonstrated that the tidal shear plays an important role in the dynamics of collapsing 
regions (Dubinski 1992; Bond & Myers 1993; Bertschinger & Jain 1994; van de Weygaert & 
Babul 1994) and therefore requires a non-spherical analysis. Non-sphericity is also necessary 
in order to explain the origin of galactic rotation through the coupling of galaxies to tidal 
torques during their collapse (Hoyle 1949; Peebles 1969; Doroshkevich 1970; Efstathiou & 
Jones 1979; White 1984; Hoffman 1986; Barnes & Efstathiou 1987; Ryden 1988; Quinn & 
Binney 1992; Warren et al. 1992). 

In this paper we construct a new analytical model to study the non-linear collapse of non- 
spherical regions in a Gaussian random field of initial density perturbations. The collapsing 
system is approximated as a triaxial ellipsoid of constant overdensity that is affected by 



2 



its own gravitational field and by the external tidal shear. The external tidal force can be 
expanded as a multipole series, with the quadrupole being the dominant relevant term for the 
internal dynamics of the ellipsoid. Since the quadrupole force and the ellipsoid self-gravity 
are linear functions of the spatial coordinates, they maintain homogeneity of the ellipsoid 
at all times. The quadrupole shear from the external density field is calculated by dividing 
the background mass distribution into spherical shells that move only radially. Under these 
approximations, we reduce the full dynamical equations of motion to a set of second-order 
ordinary differential equations. The initial conditions can be derived from a Gaussian random 
field of density perturbations with some power spectrum, and the initial peculiar velocities 
are chosen to be those of the growing mode. The equations of motion we obtain reproduce 
the linear regime behaviour of the density field. By integrating these differential equations 
we are able to describe the collapse of a triaxial sheared region into a virialized object. 
In difference from the spherical approach, the ellipsoid model can be used to study the 
influence of the external tidal shear on the triaxial collapse of overdense regions. Through 
many random realizations of the initial density field the model can examine the statistical 
properties of collapsing systems, including their triaxial shapes, orientation relative to the 
background density field, and total angular momenta, for different cosmological models. 

The outline of this work is as follows. In §2 we review the equations of motion for a 
self-gravitating ellipsoid and the approximation (Icke 1973) that allows us to extend this 
treatment to a homogeneous ellipsoidal overdensity in the universe. We also show that the 
quadrupole shear can be added to this formulation. In §3 we present our approximation for 
the non-linear evolution of this shear. The restrictions placed upon the initial conditions for 
the model are described in §4. Given a Gaussian random field of initial density perturbations, 
we derive in Appendix A the joint probability distribution for the necessary initial conditions. 
In Appendix B we describe our treatment of a collapse along a single axis, and in Appendix 
C we relate the magnitude of the external shear to the amplitude of density fluctuations on 
a given mass scale in a fashion independent of the initial power spectrum. The statistical 
properties of collapsing regions are then analyzed in §5, by applying the above model to 
many realizations of random initial conditions. Finally, in §6 we discuss the applications and 
limitations of the model and indicate potential future work. 
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2. ELLIPSOID MODEL FOR A COLLAPSING REGION 

We begin from a small density perturbation on top of a smooth background assumed 
to be a Friedmann-Robertson- Walker (FRW) universe. We assume that the background is 
composed of collisionless cold particles and allow for a no n- zero cosmo logical constant. We 
pick the origin to be near the center of a region that has a sufficiently high density so that it 
collapses before its environment. Since we are interested in the properties of the collapsing 
region, we separate the universe into two disjoint parts, the collapsing high density region 
and the rest of the universe. The boundary between these parts is taken to be a sphere 
centered at the origin. We denote the background density by p& and let the density as a 
function of space be p(x), where x is the position vector. We then define <5p(x) = p(x) — p& 
and <5(x) = 5p(x)/p^. 

The gravitational force exerted on the spherical region by the rest of the universe can be 
calculated by expanding the external gravitational potential as a multipole series (Binney 
and Tremaine 1987). Taking the sphere to have radius R, we find 

$ W = ^T[atmYt m x e , (1) 

£,m 

where x = |x|, the Yi m are spherical harmonics, G is Newton's constant, and 

ae m = Pb J d 3 sY? m 6(s)s- e -\ (2) 

|s|>i? 

The time-dependent magnitude of the external potential will be calculated in §3. For now 
we consider the effect of this potential on the inner region. 

We are primarily concerned with the acquisition of angular momentum and the shearing 
of the inner region. For these purposes, we focus on the quadrupole (I = 2) terms. The 
£ = term produces no force, and the dipole (I = 1) terms produce a uniform acceleration 
that may move the whole inner region but does not alter the shape or induce any rotation. 
This motion can indirectly affect the object because the surrounding material will have a 
new angular distribution relative to the displaced object, requiring a re-calculation of the 
multipole expansion coefficients. However, if the dipole is generated on large scales, then 
the object and its entire neighborhood move together as a bulk flow and the changes in the 
angular distribution of matter will be very small, allowing us to ignore the £ = 1 terms. 
Previous work has found the quadrupole terms to dominate the higher (£ > 3) terms (Quinn 
& Binney 1992), and so we ignore all but the £ = 2 terms. 
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We model the inner region as a homogeneous ellipsoidal overdensity (Lynden-Bell 1964; 
Lin, Mestel & Shu 1965; Zel'dovich 1965; Icke 1973; White & Silk 1979). By this, we 
mean that <5(x) is a constant inside the ellipsoid and zero outside of it. Homogeneous 
ellipsoids are advantageous because their gravitational potentials depend quadratically on 
the coordinates and therefore preserve their homogeneity and because they have non-zero 
quadrupole moments which couple to the t = 2 terms of the external potential. Obviously, 
this approximation ignores all the complexities of shapes and substructure actually present 
in the central region. However the object cannot be torqued by its sub-components, and 
the time dependence of the quadrupole moments may be well captured by the ellipsoid 
approximation. 

A convenient way to analyze the motion of homogeneous ellipsoids was discussed by 
Peebles (1980, §20). The position of the ellipsoid is given by 

r a = A a P X P, (3) 

where x is a vector inside the unit sphere and repeated indices are summed. The matrix 
A is constant in space but depends on time. The mass is evenly distributed over this unit 
sphere and thus the ellipsoid has a uniform density. If the columns of A are orthogonal, 
then they are the axes of the ellipsoid, with the lengths of the axes being the magnitude of 
the respective column. However, since in general the columns of A are not orthogonal, the 
axes may be found from the relation defining the outer shell of the ellipsoid, x a x a = 1. The 
equation of the ellipsoid is therefore 

v T A- 1T A- 1 r = 1. (4) 

To rotate to the principal axis frame we diagonalize the matrix AA T as QAQ T , where Q is 
orthogonal and A is real and diagonal because AA T is symmetric. Then we rewrite equation 
(4) as 

{Q T v) T h-\Q T v) = 1. (5) 

Thus, the axes of the ellipsoid are in the columns of Q and the corresponding axis length is 
the square root of the corresponding diagonal element of A. 

Since the gravitational potential of such an ellipsoid may be written as a quadratic 
function (Peebles 1980), we are prompted to consider a general quadratic potential 

$(r) = (6) 

where the matrix <3> a ^ is a function of the matrix A. Then the force per unit mass is just 



■5 



— V$(r), leading to the momentum equation, 



dt 2 

By substituting equation (3) in (7) we find 

d 2 A aP 
~d£~ 



<,2 '' a = (7) 



_$«7^7/3. ( 8 ) 



This is our basic equation of motion. It is crucial that the potential be quadratic in the 
coordinates, since this leads to forces that are linear in space. The magnitude of r cancels 
out, so that all the similar ellipsoidal shells behave in the same way and the ellipsoid remains 
homogeneous. 

We now must construct this quadratic potential. First, we consider the contribution 
from the ellipsoid. For an isolated ellipsoid of mass M e , the potential is given by elliptic 
functions (Peebles 1980) so that, 

$ e// (r') = \GM e [(rtfRDial a|, a 2 ) + (r' 2 ) 2 R D (al a|, a|) + (r' 3 ) 2 R D (al a|, a|)] , (9) 

where aj are the lengths of the semi-axes and r' is in the principal axis frame. The function 
Rd is defined as (Carlson 1977; Press et al. 1992), 

00 

R D (x,y,z) = -[ dt (10) 

27 ( t + z) v /(t + x)(t + y)(t + z) 

Rotating to the original frame by means of r = Qr' yields the matrix 



®ell = GM e 



Qdiag(^(A 22 ,A 33 ,A 11 ), J RD(A 11 ,A 33 ,A 22 ), J Ri 3 (A 11 ,A 22 ,A 33 ))g r . (11 



Although the potential inside a homogeneous ellipsoid is quadratic, the potential outside 
is rather complicated. This causes the smooth background outside the ellipsoid to warp, 
producing a non-quadratic potential inside the ellipsoid and causing it to become inhomo- 
geneous. In order to avoid this problem, we assume that the background remains smooth, 
with a density equal to that of the unperturbed universe (Icke 1973; White & Silk 1979). 
The mass of the ellipsoid, including the contribution from the background density in the 
volume covered by the ellipsoid, is taken to be a constant denoted by M. We then compute 
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the gravitational potential as the sum of two pieces. The first piece comes from the smooth 
background of density pb, yielding $ sp /j(r) = 2irGpbr 2 /3 and 

<, = |W^ (12) 

where / is the identity matrix. The second piece is associated with the remaining mass of 
the ellipsoid, and is given by § e u in equation (11) with the mass M e = M — pbV, where 
V = 47rdet(v / A)/3 is the volume of the ellipsoid. This approximation agrees well with N- 
body simulations (White 1993), because at early times the background is still smooth while 
at late times the background has a much lower density than the ellipsoid. 

The evolution of the background density pb oc a~ 3 in a matter-dominated universe can 
be obtained from the FRW equation for the scale factor a = 1/(1 + 2;), 

-^J = Qa^ 1 + Q A a 2 + Q R . (13) 

We define f2 = 8irGpb/3HQ, tt\ = A v /3Hq, and = 1— f2— il\ to quantify the contributions 
of non-relativistic matter, the cosmo logical constant A v , and the space curvature to the 
expansion of the universe. These quantities are evaluated at the present time, and Ho = 
d(\na)/dt\ a= i is the present Hubble constant. 

Next, we note that the £ = 2 terms of the external shear [Eq. (1)] also produce a 
quadratic potential. Manipulating the spherical harmonics in equation (1) we find 



( 2y / 6Re a22 — 2a2o — 2-\/6Ima22 — 2\/6Rea2i\ 

—2^^022 — 2v / 6Rea22 — 2<220 2v^6Ima2i 
\ — 2y / 6Rea2i 2-\/6Ima2i 4a2o / 



(14) 



Here the relation a 2 (_ m ) = (— l) m a2m was use d to eliminate m < 0. We choose the real and 
imaginary parts of 022 and 021 as well as 020, which is always real, as the five independent 
real values in the t = 2 decomposition. We will discuss the time dependence of the 02m 
coefficients in §3 and the issue of how to pick them initially in Appendix A. 

With the equation of motion (8) and the potential matrix = ($^ + + ^h ear )i 
we may evolve the system of a homogeneous ellipsoid on a smooth background undergoing a 
time-dependent quadrupole external shear as a set of nine second-order differential equations. 
For a non-zero cosmological constant, one should also add: §vac = — |A„/ Q ^. In order to 
apply this approach to a collapsing perturbation, we must relate the initial conditions of the 
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ellipsoid and the time- dependence of the external shear to properties of the initial density 
field. As described above, we center the origin of the coordinate systems on a high density 
region and construct a sphere about this point to distinguish the collapsing object from its 
environment. We consider the system at high redshift, so that (5 < 1. The coefficients 
a2m that describe the external potential as a function of the initial density field are given 
in equation (2). We pick the ellipsoid to match the average density, mass, and quadrupole 
moments of the inner spherical region at the initial time. This choice is independent of initial 
time to leading order in the linear regime. We define the average overdensity of the inner 
region at the initial time as 



m= \*n&) J dr6(r) - (15) 

\r\<R 

The mass of the region is M = (47r/3)pb-R 3 [l + S(R)], and the quadrupole moments q2 m are 
defined as 

q2m = Pt J d\5{v)r 2 Y 2 * m . (16) 

|r|<i? 

To match these quantities, we consider an ellipsoid with semi-axes ci, c 2 , and C3 oriented at 
some angle relative to the coordinate axes. We pick the overdensity of the ellipsoid to be 
S(R). Then the mass of the ellipsoid is (47r/3)p&[l + 8{R)\c\C2C^, which when matched to 
the mass M of the inner region gives 

C1C2C3 = R?. (17) 

Finally, the quadrupole moments of the ellipsoidal overdensity are chosen to match those of 
the actual inner region. First, we label the points of the ellipsoid by r = QCx, where x is a 
point inside the unit sphere, C = diag(ci, C2, C3), and Q is an orthogonal matrix whose jth 
column is the direction of the axis with length Cj. We then define the matrix N = QC 2 Q T 
and note the following integral over the volume of the ellipsoid, 

Pb$ J d 3 rnrj = p & 5cic 2 c 3 J d 3 x(Q irn C mm x m )(Q jn C nn x n ) 



ellipsoid I X I<1 



, rp3 l ^5 mn \ -MlN-- 



(18) 



where M e = (4n/3)p b 5( y R)R 3 is the mass of the ellipsoid above the background. Now we 
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may find the q 2 m of the ellipsoid in terms of the matrix N; for example, 



Req 22 =p b 5 J ^V^( r i- r 2 2 ) = M Vl6L (Aril_iV22) - (19) 

ellipsoid 

These relations may be inverted, so that given all the q 2m we find N to be 

/ Re q 2 2 ~ ^<?20 -Im ?22 -Re q 2 \ \ 

-Img 2 2 -Reg22--^g20 Img 2 l 



/40tt 1 
" 3 M P 



+ tI, (20) 



\ -Reg 2 i Im<72i 4?<?20 / 



where the constant r in front of the identity matrix is unknown. We then diagonalize N 
to find Q and C 2 , the latter depending on r, and impose the condition in equation (17) 
to determine r. With this condition we find the lengths and directions of the axes of the 
ellipsoid. We then set the initial value of A to be QC. 

Since the equation of motion is second-order, one must specify the initial velocities. We 
pick the velocities so that the density field is a pure growing mode, consistent with the fact 
that we normalize the power spectrum today when only the growing mode had survived. In 
linear theory the peculiar velocities are given by (Peebles 1980, §14) 

5v = amy (21) 

where g is the peculiar gravitational acceleration, H is Hubble's constant, and Q = pb/pc 
is the cosmological density parameter. Here, / = (a/ D)(dD / da) as f2 - 6 , a is the expansion 
factor of the universe, and D is the linear growth factor. Since g = — $r = — $Ax, the initial 
velocity field is 

d A If 

v =M*= Hr -im* r ' (22 > 



and so we find 



^— = HA a/3 - ^J-^A^. (23) 
dt 3HQ v ; 



In this way we have related the initial conditions of the equation of motion to 5 and the five 

After evolving the ellipsoid through the matrix A, we wish to measure properties of 
the ellipsoid, such as its angular momentum and energy. The axes of the ellipsoid are the 
eigenvectors of AA T and the lengths of the semi-axes are the square roots of the eigenvalues. 
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Labeling these lengths as Cj, the potential energy of the ellipsoid is given by (Binney and 
Tremaine 1987), 



oo 

W = --GM 2 
5 



oo 

/ . ^ (24) 



Here we have neglected the potential energy from the background density and from the tidal 
shear, both of which are small at late times when this quantity is of interest and when 
M e « M. The kinetic energy equals (Peebles 1980, §20), 

1 f 2,s PedetA f dA ik dAi k , , M fdAdA T 



T= -p e \ v z d 6 r= ^ " / , xV = — tr , (25) 

2^ e / 2 / dt 10 \dtdt ]l 



|x|<l 



where p e is the ellipsoid density. The total energy is = T + W. Next, the jth component 
of the angular momentum is 

/r J Amp jyj d A mn 

(r x v)jd 3 r = p e det A / e ^ km A kn ^—x n x p = ™ e 3km'?fL_A kn , (26) 

|x|<l 

where e> km is the Levi-Civita anti-symmetric tensor. The mass of the ellipsoid M is con- 
stant, and at late times the contribution from the smooth background is small. With these 
definitions we can construct the spin parameter A that measures the amount of rigid rotation 
acquired by the ellipsoid before virialization, 



A ee (27) 
GM 5 / 2 v ; 



3. TIME DEPENDENCE OF THE EXTERNAL SHEAR 

In the equation of motion (8) for the ellipsoid, one must specify the time dependence 
of the external shear. The model described in §2 considers only quadrupole tidal forces 
and therefore requires a limited amount of information about the mass distribution outside 
the ellipsoid boundary. In this section, we will try to develop a simple model for the time 
dependence of the a2m coefficients. 
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The simplest approach is to treat the entire volume outside the spherical boundary of the 
inner region by linear perturbation theory. The motivation for this approach is that because 
the inner region is a high density peak, it is likely to be the first object in its neighborhood to 
collapse. Other high peaks are sufficiently far away that their collapse does not substantially 
change the tidal force, since such collapses do not move matter on large angular scales as 
seen from the origin. If the tidal field is determined by scales much larger than the inner 
region, the evolution of the quadrupole moments would be well treated by linear theory. 

In linear theory, the fractional overdensity <5(x) scales by a uniform growth factor D(t) 
when considered in comoving coordinates (Peebles 1980). /,From equation (2), we have 



where the s coordinates are in physical (non-expanding) space. Rescaling to comoving co- 
ordinates, however, makes no difference because the magnitude of s appears equally in the 
numerator and denominator. The density pb scales as (1 + z) 3 , where z is redshift. Thus in 
linear theory, 



where t{ is the time when we pick the initial conditions. In a flat matter-dominated universe, 
D oc (1 + z)^ 1 and a2 m oc (1 + z) 2 . 

The time dependence given in equation (29) can be substituted into the equation of 
motion for the ellipsoid. When this is done and an ensemble of random ellipsoid realizations 
are integrated until full collapse, we find the average value of the spin parameter A to be 
smaller than 0.01. Since N-body simulations typically find (A) ss 0.05 (Barnes & Efstathiou 
1987; Warren et al. 1992), the linear approach apparently underestimates the external 
torque. It is not surprising that the external shear cannot be purely described by linear 
theory, since the material just outside the boundary of the object has a density similar to 
that of the object. As the ellipsoid collapses, the nearest shell surrounding it follows its infall. 
Since the contribution of a mass element to aim depends on the inverse cube of its distance 
from the origin, the fact that material near the object tends to be closer than its initial 
comoving position means that it should produce more torque than linear theory predicts. 

Based on the fact that the material around the peak tends to be closer than linear 
theory predicts, we divide the external mass distribution into N spherical shells, all centered 




(28) 



|s|xR 




(29) 
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on the origin, and evolve each shell according to the mass interior to it (Gunn & Gott 1972). 
Each shell carries a fraction of the quadrupole shear, initially determined from the Gaussian 
random field, and this shear is scaled according to the radial collapse of the shell. The 
essential approximation is that we neglect the tangential redistribution of the matter within 
shells. This approximation is justified so long as the motion of the matter subtends a small 
angle, as viewed from the origin, which is always true for matter at large radius. Using 
growing mode velocities, we need to specify only the average overdensity interior to the shell 
S in order to determine its motion. 

To calculate the time dependence of the shear in this approximation, we first note that 
a 2m depends on <5p(x) but not on rescaling of the radial coordinate. We find Sp by solving 
the spherical evolution from some initial overdensity S(ti) using the FRW equation (13), 

d 2 r _ 4:7i G pi A v r 

with initial density pi — (1 + S)pi,(ti) and r(ti) = 1. The initial velocity is to leading order 
H(ti)(l-5/3), with 

H(ti) = H ^n(i + zif + n R (i + zif + n A , (31) 

where Zi is the redshift at the initial time ti. Thus for a given S, we may integrate equation 
(30) to find r(t). We then find a2 m {t) oc 5p by considering the difference between the top-hat 
density, p sp h oc r~ 3 , and the background density, p& oc (1 + z) s . Thus, 

°2m(*) = (l + 6i)r(t)-*-([l+z(t)]/[l + Zi])^ 

02m(ti) Si 

For a sufficiently early initial time in the matter-dominated epoch, Si <^ 1 and Q ~ 1, so 
that equation (32) matches the linear theory result of aim oc (1 + z) 2 . 

To implement this approach, we wish to consider iV shells, with radial boundaries R a , 
for a — 0, 1, . . . , N. We take Rq to be the radius of the spherical boundary that separates the 
ellipsoid region from the external region and Rjy = oo to indicate that the outer shell includes 
everything outside the outermost spherical shell. Initially each shell carries a fraction of the 
quadrupole shear 

a£2 = Pb J d 3 sY 2 * m S(s)s-\ (33) 

i?„-i<|sj<i?„ 

where n = 1,2, ... ,N labels the shell. We pick the characteristic overdensity of each shell 
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as 

S {n) ^^(6[R n ^]+5[R n ]). (34) 

—<n) - 

For thin shells the differences between the definitions 5 and 5(R n ) defined in equation 
(15) is small. The overdensity surrounding the outermost shell Rn-i is negligible, so we set 

-(AO 

5 =0 and treat the tide of the outermost shell using linear theory. 

The above model for the exterior region completes the necessary set of equations and 
allows us to evolve the ellipsoid beyond the linear regime according to equation (8). This 
approach requires a large set of initial data: the overdensities 5 for all the shells and for the 
inner region, the quadrupole shears d^ m for all outer shells, and the quadrupole moments qi m 
for the inner region. This set comprises 6N + 5 real numbers, randomly drawn in a highly 
correlated way from a Gaussian random field of initial density perturbations. By treating 
the random density field in spherical coordinates, as suggested by Binney and Quinn (1991), 
we can separate this set of random numbers into six independent sets, each of which is 
distributed as a multi-dimensional Gaussian distribution. We describe how to find these 
distributions in Appendix A. We may then draw the initial conditions from the proper 
probability distribution by taking linear combinations of 6N + 5 random Gaussian deviates. 



4. COLLAPSE AND CONSTRAINTS 

If the ellipsoid is simply evolved forward in time, one soon reaches the obvious problem 
that the axes collapse at different times. Typically, one axis turns shorter and collapses 
first, forming a pancake (Lin, Mestel, Shu 1965; Zel'dovich 1970; Peebles 1980) in which the 
baryons shock and the dark matter goes through violent relaxation. The computation of 
angular momentum should not end at the time of the short axis collapse, since the other 
two axes are still extended and give the object a significant quadrupole moment with which 
it couples to the tidal shear. Previous work on spherical objects (Ryden 1988) has shown 
that most of the angular momentum is acquired near turnaround, since most of the time is 
spent at this phase. In the ellipsoid model, turnaround is not a single epoch and it appears 
necessary to follow the system as long as the quadrupole moments are large, i.e. until the 
long axis turns around. After the short axis collapses, it makes a small contribution to 
the quadrupole moment of the ellipsoid. The quadrupoles are proportional to the difference 
between the squares of the axis lengths and so the exact dynamics of the virialization process 
along the short axis has little effect on the acquisition of angular momentum. 
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We therefore follow the approximation of Bond and Myers (1994) and impose a strict 
cutoff on the collapse of an axis, namely that no axis may collapse below 40% of its maxi- 
mum length. This keeps the dynamics from approaching the singularity at zero length and 
simulates the virialization of the corresponding axis. The 40% cutoff was picked to allow the 
ellipsoid to be slightly flatter than the spherical virial theorem would predict, but because 
the collapsed axis is much shorter than the non-collapsed ones, the exact cutoff value does 
not affect the angular momentum acquisition significantly. The implementation of this con- 
straint is non-trivial because the ellipsoid may be rotating and the columns of the matrix A 
are in general not orthogonal. The details of this implementation are presented in Appendix 
B. Note that the vorticity of the ellipsoid is conserved under the equation of motion (8), but 
is not conserved through this treatment of the collapse of an axis. Also, since the kinetic 
energy of the ellipsoid is being altered, energy is not conserved. Because the change in the 
total energy during the dynamics is usually small, we save the last value of the energy before 
the first axis collapses and use it in the calculation of the spin parameter A. 

In principle, we would like to evolve the ellipsoid until all three axes collapse, by which 
time the angular momentum gained by tidal torques (i.e. not by accretion) will be complete. 
The first difficulty is that in some cases one or two axes do not collapse at all. The existence 
a high-density region does not guarantee that it will not be sheared by its environment to 
form a wide pancake. This would be true even if we were to use linear theory for the shear, 
as can be seen from the Zel'dovich approximation (Zel'dovich 1970). In this approximation, 
if a volume element begins with an outward peculiar velocity along one of its axes then it 
will continue to move outwards for all later times. Because the region is overdense, all the 
growing mode velocities resulting from the ellipsoidal perturbation will be inward. But since 
the trace of Q s hear is zero at least one eigenvalue of & s hear is negative, and therefore in at 
least one direction, the peculiar velocity will be less inward than for the isolated ellipsoid. We 
wish to study regions that end as spatially bounded objects. In f2 = 1 cosmologies, following 
the Zel'dovich approximation, we require that the radial component of the initial peculiar 
velocities inward in all directions. This corresponds to requiring that the potential matrix 
be positive definite. In Q < 1 cosmologies, however, this condition is not sufficient because 
mildly overdense regions may still expand forever. Only regions which exceed the critical 
density collapse, so we instead require that the inward radial component of the peculiar 
velocities exceed that of the growing mode of a spherical top-hat perturbation at the critical 
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density, which corresponds to an initial overdensity 
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Due to nonlinear effects of the external shear, our requirement is not a perfect discriminator 
as to whether an axis will expand forever, but it at least matches the weakly non-linear 
theory. The consequences of this constraint are examined analytically in Appendix C. 

Since we base the time dependence of the shear on the density of the spherical shell 
collapse, the shear force diverges as the shell collapses to a zero radius. Because the average 
overdensity of the innermost shell is not so different from that of the ellipsoid, the two objects 
collapse at similar times. But in most cases, the action of the shear causes the long axis 
to collapse significantly later than the inner external shell. We must end the integration at 
some time before the collapse of the inner tidal shell. We therefore choose to end at the time 
at which a sphere whose overdensity was that of the initial ellipsoid would collapse to a zero 
radius. We then constrain the initial conditions so that none of the exterior shells has an 

— in) 

overdensity 5 greater than 95% of the initial density of the ellipsoid. This ensures that 
the external tidal shells do not collapse before the integration ends. This condition has the 
important property of tending to make the inner region a density peak; we do not explicitly 
require that it is truly a peak, but if the overdensity as a function of radius is forced to drop 
as one moves away from the central region, only extreme variations as a function of angle 
would allow the density to rise radially in a particular direction. 

The last and most significant constraint we impose is that the inner region must have a 
high overdensity. This condition underlines the separation between the collapsing region and 
the rest of the universe. We implement this constraint by requiring S(Ro) > VminC, where 
v-min = const and er is the r.m.s. amplitude of mass fluctuations 5M/M for such regions. As 
usual, the latter is calculated to be (Bardeen et al. 1986), 



where P(k) is the power spectrum and ji is a spherical Bessel function. In §5 we will typically 
pick Vrnin ^ 2, thus forcing the central region to be rare. 

The above three constraints — (i) 5(Rq) > v m i n o\ (ii) the overdensity of the surrounding 
shells being less than 95% of 5(Rq); and (iii) all peculiar velocities being inwards — are not 
implemented directly into the probability distributions for the initial conditions as derived 
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in Appendix A. Rather, we generate sets of the initial conditions and then reject cases that 
do not satisfy these constraints. Since the strictest condition is S(Rq) > u m i n a, we set up the 
probability distribution so that S(Rq) depends on only one Gaussian deviate, which allows 
us to reject a set of initial conditions on the basis of one random number. Then we can 
choose the overdensity of the innermost shell to depend on only one more deviate, allowing 
another straightforward test. 

We call the fraction of sets of initial data that meet the above conditions the acceptance 
rate A. Prior to imposing these constraints, the Gaussian random field was unconstrained. 
Because the central region is at a mass scale M = (47r/3)i?Q, the number density of such 
regions is hq = po/M, where po is the mass density today of non-relativistic matter (i.e. 
the matter making up the ellipsoid). This means that after applying the constraints, the 
number density of accepted regions is Auq. If we now construct an ensemble of such regions 
and, after evolving each of them with our model, find that a given property occurs in some 
fraction / of them, then the actual number density of such objects is predicted to be f An$. 
We apply this approach elsewhere to predict the number density of black holes form in the 
initial collapse of objects with very low values of angular momentum (Eisenstein and Loeb 
1994, hereafter Paper II). 

This completes the description of our model. We next study the statistical properties 
of collapsing regions in specific cosmological models, defined by the values of the cosmolog- 
ical parameters: Q, A v , and Hq, as well as by the power spectrum of primordial density 
perturbations. For a given mass scale, we consider many realizations of the random initial 
conditions, evolve each realization separately, and then combine the results to analyse the 
statistical distribution of shapes and angular momenta of collapsing regions. 



5. RESULTS 

To illustrate the behaviour of the model, we adopt a cold dark matter (CDM) power 
spectrum (Bardeen et al. 1986) of initial density perturbations for two different universes: 
a flat universe with Q = 1, and an open universe with Q = 0.2; both cases assuming A v = 
and Hq = 50km/s/Mpc. We normalize the power spectrum by choosing the present r.m.s. 
amplitude of mass fluctuations within an 8/t _1 Mpc radius sphere to be 1/6, where b is the 
bias parameter. We pick the initial time to correspond to redshift z% = 1000, and then scale 
the power spectrum using the growth factor D(to)/D(ti). 
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The parameters of the model are the mass of the central region M, the minimum peak 
height Vmini and the radii of the external shells. We consider a variety of mass scales between 
10 8 and 10 15 M . As peak thresholds we use both v m i n = 2.5 and 2.0. The radius of the 
innermost spherical boundary Ro is fixed to match the mass scale being studied, and we pick 
all other radii to be multiples of this inner radius. We choose twenty shells, with boundaries 
at 1.5, 1.75, 2, 2.5, 3, 3.5, 4, 4.5, 5, 6, 7, 8, 9, 10, 12, 15, 17, 20, and 30 times the inner radius. 
The only important number in this list is the first one, which implied that the innermost 
shell is thick. This is done in order to keep the acceptance rate reasonable, as the overdensity 
of the closest outer shell must be less than 95% of that of the inner region. Were this shell 
to be thin, the overdensity constraint would only rarely be satisfied. The finite thickness of 
this shell ultimately limits the radial resolution of the external shear profile. 

We consider the flat universe with b = 1 for a variety of mass scales and for v m i n — 2.5. 
The average value of A increases slightly with mass, from 0.034 at 10 8 M to 0.038 at 10 12 M 
and to 0.049 at 1O 15 M . Figure 1 shows the probability distribution of the spin parameter 
-P(A) for the mass scales of 10 12 M and 10 15 M in a flat universe. The histograms are 
based on samples of random realizations that result in 2 x 10 5 accepted systems for each 
mass scale. The weak dependence of the mean values of (A) on the mass scale is in good 
agreement with with N-body simulations (Barnes & Efstathiou 1987; Warren et al. 1992). 

The primary reason for the actual dependence of (A) on M is the shape of the power 
spectrum. With the cold dark matter power spectrum, higher mass scales have larger average 
quadrupole moments (relative to a MR 2 ). This results in more elongated initial ellipsoids 
and therefore in stronger couplings to the external torque. In addition, the fixed choice 
of shell boundary radii introduces a systematic mass dependence because the higher mass 
systems have a steeper density profile and a larger fraction of the shear carried by the inner 
shells (these two profiles are related, as shown in Appendix C). Consequently, the effective 
shell resolution appears coarser for high mass objects. The limited resolution results in a 
slight underestimate of the shear at late times, when the shear amplitude is large. 

For Vmin = 2.5 only a fraction 0.00621 of all regions satisfy the constraint S(Rq) > fmin^- 
The remaining two constraints (shell densities and inward velocities) are satisfied 64.2% of 
the time for 1O 8 M and 89.7% for 1O 12 M . Most of these additional rejections are due to 
the inner shell having too high an overdensity; rejections due to the peculiar velocities not 
being inward occur in less than 3% of the cases (cf. Appendix C). With peak threshold values 
of 

v-min — 2.5 and 2.0, the average values of v are 2.84 and 2.40 respectively, independent of 
mass scale. 
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For a flat universe, the model is strictly independent of the bias parameter or the Hubble 
constant. Both parameters rescale space and time without altering the collapse dynamics. 
This does not hold in an open universe, where an intrinsic time scale appears at the transition 
point to an open curvature-dominated expansion of the background. 

The value of A is found to be anti-correlated with v — 8/a, so that higher peaks tend 
to have lower angular momentum. The linear correlation between the two is about 0.21. 
Therefore, the choice of v m i n affects the distribution of A. For a flat universe and M = 
10 12 M we find (A) = 0.038 for u min = 2.5 and (A) = 0.051 for u min = 2.0. For 10 15 M 
we find (A) = 0.049 for v m i n = 2.5 and (A) = 0.069 for v m i n = 2.0. If we use a particular v 
for all iterations, as opposed to the usual requirement of v > v m i n , then for 10 15 M we find 
(A) = 0.089 for v = 2, (A) = 0.060 for v = 2.5, and (A) = 0.042 for v = 3. 

In all of the accepted objects, the shortest axis reaches collapse by the end of the integra- 
tion (i.e. by the time a spherical perturbation with the same initial density as the ellipsoid 
would have reached zero radius). For v > v m i n = 2.5, the middle axis reaches collapse in 
78% of all iterations, while the long axis reaches turn-around in 84% of the iterations and 
reaches full collapse in 1.5% of all cases. If we force v — 2, the corresponding numbers are 
65%, 58%, and 0.6% respectively for M = 1O 12 M , with a weak dependence on mass. The 
failure of the long axis to collapse is expected; since it is being pulled out by the shear, it col- 
lapses later than the unperturbed spherical case. The few times that the long axis happens 
to collapse result from the small window of opportunity between the sphere reaching 40% 
of its turn-around radius (defined as axis collapse) and reaching a zero radius. The failure 
to turn-around usually is the result of the shear leaving one axis with only a tiny inward 
velocity so that it undergoes a very long excursion, but it can also result from extreme cases 
where the rapidly increasing shear toward the end of the integration causes one axis to ex- 
pand again after it had already been contracting. The dependence of the amount of collapse 
on v reflects the tendency of the shear to pull lower peaks apart. Although we require that 
the velocities are initially inwards, this may not be sufficient to guarantee that the object 
remains bounded under the influence of the nonlinear shear. As suggested in Appendix C, 
the constraint to have inward velocities is almost negligible for v > 2.5, with acceptances 
over 97%, but is important for v — 2, where the acceptance drops to 82%. The drop means 
that more systems have one axis with a nearly zero peculiar velocity; these axes are very slow 
to turn-around. Low-u systems therefore have larger axis ratios resulting in larger couplings 
to the quadrupole torques. The correlation between the peak height v and the axis ratios 
leads to the above dependence of A on v. 
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Let us now turn to discuss the axis ratios of the evolving ellipsoids. Figures 2-8 present 
the axis ratios b/a and c/a for various ensembles of ellipsoids. The semi-axis lengths a, b, 
and c are ordered so that a > b > c. Panels (a)-(e) show the ensemble at different times, 
while panel (f ) shows the mean of the points in the previous panels. 

Figure 2 presents the evolution of the axis ratios for a random ensemble of 10 15 M 
ellipsoids with v = 2 in the Q = 1 cosmology. We chose to fix v rather than use a threshold 
value Vmin so that the integration of each ellipsoid would end at the same time. The time 
slices in the plots show each ellipsoid at a constant fraction of the collapse time t max for a 
spherical top-hat perturbation of the same initial density. We track 5000 ellipsoids and show 
the axis ratios at the initial time and 25%, 50%, 75%, and 100% of t max - We set b = 1.3 so 
that the end time is close to z — 0. 

The most significant aspect of this plot is the progression of the points from a quasi- 
spherical initial state to a prolate end state. At early times [panels (a)-(c)] the short axis 
falls behind in the expansion so that the points move down the plot. At the time of panel (d) 
the short axis has generally collapsed; however the typical resulting shape is not an oblate 
pancake or a prolate filament but rather triaxial. As the evolution continues the middle axis 
turns around and collapses, so that most objects become prolate by the end time in panel (e). 
Note the slight regeneration of high c/a points between panels (d) and (e); this corresponds 
to the beginning of a collapse along the long axis. Indeed, 37 of these 5000 ellipsoids reached 
long- axis collapse (i.e. the long axis reached 40% of its maximum length). 

In contrast to previous discussions on collapsing ellipsoids (Lin, Mestel & Shu 1965; 
Zel'dovich 1965; Icke 1973; White & Silk 1979), we find that the collapsing region geometry is 
primarily determined by the external shear and not by the initial anisotropy of the ellipsoids. 
To demonstrate this result, we alter the model and replace the initial ellipsoid by a sphere. 
The results for a mass of 1O 15 M , Q = 1, and v — 2, are shown in Figure 3. The evolution is 
very similar to Figure 2 despite the qualitative change in the initial conditions, although the 
ratio of the long axis to the short axis is not quite as large at late times. Furthermore, even 
in the full ellipsoid model, the direction of the long axis is dominated by the direction of the 
shear rather than the direction of the long axis of the initial ellipsoid. This is particularly 
the case for low mass scales and for higher u, where the initial quadrupoles are weak and the 
initial conditions are closer to the spherical example. Even an initially spherical object can 
gain angular momentum because the shear direction changes as the relative weighting of the 
shells evolves with time. The resulting value of A is lower for the spherical case, however, 
indicating that the anisotropy of the initial ellipsoid is important for the acquisition of 
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angular momentum. For the parameters in this paragraph, (A) = 0.089 when the initial 
object is the usual ellipsoid, but (A) = 0.028 when the initial object is changed to a sphere. 

In Figure 4, we consider the mass scale of 10 12 M and v = 2 in the Q = 1 cosmology (the 
redshifts assume 6 = 1.3). Compared to Figure 2, the initial ellipsoids are more spherical 
and have weaker quadrupole moments. However, the axis ratios at late times are similar to 
those at the 10 15 M mass scale. 

We next consider the open cosmology with = 0.2 and 6 = 1. For the 10 15 Mq mass 
scale we pick v = 2.2 so that the final redshift of the calculation matches that of the flat 
universe case. The results are presented in Figure 5. Here we have chosen the time slices 
so that the redshifts are the same as in Figure 2. The slices are no longer evenly spaced in 
time because of the different cosmology. Due to changes in the power spectrum, the objects 
are initially more spherical in an open universe, but they end with very similar axis ratios 
to those obtained in a flat cosmology. 

A major difference between the open and the flat cases is that the constraint on the 
initial velocities is far more severe in the open cosmology. In the open case we require that 
the peculiar velocities not only be inward but that they be larger in magnitude than those 
induced by a critically bound sphere. For the above open cosmology and the above mass 
scale this condition [cf. Eq. (35)] requires that a spherical top-hat perturbation have v > 1.1. 
Including the disruptive effects of shear (cf. Appendix C), we find that the acceptance rate 
for this constraint drops down to 23% from its value of 83% for Q = 1. Thus, most v — 2 
peaks have one axis that is not critically bound. In Figure 6, we relax the constraint on the 
velocities entirely so as to investigate how the general peak behaves in this open universe. 
There are two notable differences between this plot and either Figure 5 or Figure 2. First, the 
collapse of the shortest shown by the accumulation of the points near the horizontal 

axis, occurs in panel (c) rather than in panel (d). The relaxation of the velocity constraint 
allows larger shears to enter. This not only tends to pull out the long axis but also tends to 
push in the short axis, since the shear is a traceless matrix. The short axis in this sample 
therefore collapses faster. Second, at the final time in panel (e), the objects are significantly 
more filamentary (note that small deviations in the plotted ratio c/a near zero correspond 
to large changes in the filament shape a/c). This effect occurs because the long axis expands 
without bound and the short axis collapses earlier than before. Comparing Figures 2 and 6, 
we see that when viewed at equal redshifts, open universes have significantly more filamentary 
structure than flat universes. This occurs primarily because of the stronger influence of shear 
on the collapsing regions. 
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Next, we consider the low-spin tail of the population of collapsing regions (cf. Fig. 
1). Because of its low angular momentum, the gaseous component of the systems in this 
tail can collapse to small radii and form compact massive objects. The low-spin systems 
therefore provide environments that may favor the formation of seeds for quasar black holes; 
we discuss this possibility in detail in Paper II. For the present discussion, let us examine 
the qualitative properties of low-spin systems selected out of a population of collapsing 
regions with Vmin = 2.5 on the mass scale of 10 8 M in a flat universe (including the velocity 
constraint as usual). Figure 7 shows the full distribution of the ensemble of runs for this mass 
scale. The different panels show the ellipsoid axis ratios at fixed fractions of the collapse 
time t max of a spherical perturbation with the same initial density as the ellipsoid. Since v 
is not constant, the collapse time is different for different ellipsoids. Thus, the panels in the 
figure are not equal time slices, but instead show each ellipsoid relative to a measure of time 
scale for its development. 

We now compare this full distribution to the distribution of the subset of low-spin sys- 
tems. For the definition of a low-spin system, we use the criterion developed in Paper II 
to determine whether the gas in a system can settle into a sufficiently compact disk with 
a viscous timescale shorter than the star formation time ~ 10 7 yrs (i.e. y < 1.1 x 1CT 3 in 
the notation of Paper II). This condition translates to a maximum value for the angular 
momentum per unit mass of the system, J/M = 8.6 x 10 24 cm 2 s _1 . Since this threshold in- 
volves only the total angular momentum, it does not precisely correspond to a bound on the 
spin-parameter A (which depends on the energy as well); nevertheless, it is approximately 
the condition A ^ 1CP 3 . In a set of 2 x 10 5 accepted systems, 232 satisfy this limit. In 
Figure 8, we plot the axis ratio histories of these objects. Evidently, these objects evolve in a 
much more spherical way than the typical cases shown in Figure 7; over half of the low-spin 
systems have axis ratios between 2:1 and 1:1 at the end time. This striking result is a con- 
sequence of two correlated effects. First, the low-spin systems have closer to spherical initial 
conditions, but second and more important, the quadrupole shear in their neighborhood 
is weak. One should also consider the existence of a background population of accidental 
low-spin objects. Because of the time dependence of the shear, it is possible for the angular 
momentum of a filamentary object to be rapidly changing at late times and by chance be 
small at the end of the integration. Such elongated objects are obviously not environments 
that would favor black hole formation. While some of these background objects do occur, 
Figure 8 demonstrates that most of the low-spin events are indeed close to spherical during 
their evolution. 
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We close this section by showing four individual cases for the collapse of a 10 5 M 
object in an Q = 1 universe with b = 1.3. Figures 9a and 9b show two objects with v — 2, 
while Figures 9c and 9d show two objects with v — 3. In all cases, the first panel shows 
the semi-axis lengths as a function of time (note the different time scales between the two 
values of v). The freezing of the axis length after its collapse (cf. Appendix B) is evident. 
The second panel shows the evolution of the overdensity 5p/p = 5 for the ellipsoid (solid 
line), as well as for the spherical perturbation with the same initial density (dotted) and for 
the Zel'dovich approximation when applied to the initial quadratic potential (dashed). The 
slope discontinuities in the ellipsoid density are a result of halting the collapse of the axes; 
if this treatment was avoided the ellipsoid would have reached an infinite density slightly 
before its spherical counterpart. The Zel'dovich approximation does well at early times but 
underestimates the density after turn-around. The third panel shows the spin parameter as 
a function of time. While A grows roughly linearly with time, there are significant secondary 
variations which are more pronounced at lower mass scales. The four different figures are 
shown to illustrate a variety of collapse conditions. Figure 9a ends with relatively large values 
of A and the axis ratios. Figure 9b was picked because of the unusual wiggle in its evolution; 
such wiggles result from a time-varying shear and are not uncommon in the v — 2 systems 
when the overdensity is not high enough to dominate the shear. Figure 9c is a typical high-z/ 
object; its A and axis ratios are close to the mean. Finally, Figure 9d has a relatively low A 
but a large ratio of its long to short axes. 



6. DISCUSSION AND CONCLUSIONS 

In this work we have developed a model for the non-linear collapse of a triaxial overdense 
region out of a Gaussian random field of primordial density perturbations. The model ap- 
proximates the collapsing region as a homogeneous ellipsoid. We assume that the collapsing 

2 

mass originates in a spherical volume around a high density peak and select an ellipsoid that 
matches the mass, mean overdensity, and quadrupole moment of the initial overdensity field 
in this volume. This choice is independent of redshift to leading order in linear perturbation 
theory. The mass distribution outside the sphere exerts a tidal torque on the ellipsoid and 
spins it up. We calculate the quadrupole moment of the external shear as a function of time 
by dividing the background density field into thin spherical shells that move only radially 

2 As shown in Appendix C, only high density peaks are likely to survive as bound systems under the 
influence of the external tidal shear from their environment. This effect provides another reason for 
associating virialized objects, like galaxies or clusters, with high density peaks. 
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according to the mass interior to them. The dynamics of the ellipsoid is determined by its 
self-gravity and the external shear through a set of nine ordinary differential equations. Both 
forces are linear in the coordinates and therefore maintain homogeneity of the ellipsoid at all 
times. In Appendix A we have developed the formalism necessary to randomly determine the 
initial conditions for this model in the appropriate correlated way from a Gaussian random 
field of initial density perturbations. 

The above model was applied to a restricted set of initial conditions that are more suitable 
to its assumptions. In particular, we studied the statistical properties of rare high density 
peaks with a mean overdensity 5 2a. In a bottom-up hierarchy of structure formation, 
most objects evolve from a quasi-spherical initial state to a pancake or a filament and then 
to complete virialization. As demonstrated by Figure 3 (where the initial conditions are 
spherical), this evolution history of shapes is primarily induced by the tidal shear and not 
by the initial triaxiality of the ellipsoids. Thus, the existence of sheets and filaments in the 
universe is not a result of the Lin-Mestel-Shu (1965) instability, but rather an environmental 
effect; namely the ellipsoids are being sheared by nearby mass concentrations. As shown in 
Figures 2 and 6, the redshift evolution of the triaxiality of systems on a given mass scale can 
be used to discriminate between an open and a flat universe. However, the average value of 
the spin parameter (cf. Fig. 1) (A) 0.04, is found to be only weakly dependent on the 
object mass or the cosmological parameters, in agreement with N-body simulations. There 
is a modest dependence of A on the peak height v. 

The ellipsoid model incorporates significant qualitative improvements over previous an- 
alytical investigations of A. Other models (Ryden 1988; Quinn & Binney 1992) considered 
spherical dynamics for the collapsing region and assumed zero initial peculiar velocities rather 
than growing mode velocities. The latter assumption causes a significant overestimate of A 
since it makes the object expand to a larger radius. Quinn & Binney (1992) considered the 
dipole term of the external potential to be most important, but they worked in the instan- 
taneous rest frame of the center of mass rather than in the actual accelerating rest frame of 
the center of mass. This led to a nonzero torque even in a uniform gravitational field. 

In assessing the applicability of our model, we return to the initial ellipsoid. The matter 
inside this ellipsoid makes up the final object. The ellipsoid is picked to match the mean 
overdensity and quadrupole moment of the inner spherical region. As it traces the density, 
one may consider the boundary of the ellipsoid as an approximation to the smoothed isoden- 
sity contour of the inner region. The model then follows the evolution of this contour and 
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would therefore be useful for applications where the final object is indeed associated with 
the initial high-density contours. 

A potential example for such an application is the collapse of clusters of galaxies. Clusters 
correspond to high-density peaks on the 10 14 " 15 M mass scale. In the bottom-up hierarchy 
of structure formation the galaxies that inhabit the cluster form before the cluster collapses. 
Since galaxies tend to form in high-density regions, the high-density environment of the 
cluster peak favors galaxy formation. Galaxies therefore tend to trace the smoothed density 
field near the top of the peak, and the ellipsoid model could describe their motion during the 
formation of the cluster. Obviously the model cannot describe the final virialization process 
of clusters, where violent relaxation erases in part the signature of the initial conditions. The 
ellipsoid model predicts however that the outer unvirialized parts of clusters which are still 
infalling today would be highly triaxial with the axes ratio distribution as shown in Figures 
2 and 5. This prediction can be tested by observations of the galaxy distribution around 
clusters or by deep x-ray imaging of their surrounding gas distribution. In fact, the Virgo 
cluster is observed to be elongated considerably along the line of sight based on Tully-Fisher 
distances (Fukugita, Okamura, & Yasuda 1993). The existence of prolate systems allows for 
a systematic contamination of optical samples of rich clusters by cases of a chance alignment 
between a prolate filament and the line of sight. Similarly, clusters that are elongated 
perpendicular to the line of sight may not have been identified as clusters when examined 
by spherical filters (e.g. Abell 1958). In addition, prolate clusters introduce a systematic 
bias into the Sunyaev-Zel'dovich (SZ) effect by favoring low values of the Hubble constant 
when the data analysis is done using a spherical model for the cluster. Indeed, attempts 
to determine the Hubble constant from SZ measurements tend to get relatively low values 
of Hq, occasionally below the permitted lower bound of 50 km s _1 Mpc _1 (McHardy et al. 
1990; Birkinshaw & Hughes 1994, and references therein). Conversely, when the value of the 
Hubble constant is eventually determined by other techniques, it would be possible to get 
constraints on the triaxiality of clusters from their SZ effect. 

In some other applications, the final object evolves from a region that has little to do 
with the density contours near the top of a high density peak. For example, the matter that 
forms a galactic halo may come from regions of low initial overdensity near the peak. These 
low density regions can be assembled into the halo not only by the gravitational attraction 
of the peak itself but also by the external shear. Initially, the total densities of the "low" and 
"high" density regions are similar because the background density dominates in both cases. 
The actual boundary of a virialized object today is an equal collapse-time surface that maps 
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into a complicated surface in the initial density field according to the intervening action of the 
external shear. To illustrate this complicated situation, consider a high overdensity region 
acted upon by a quadrupole tidal shear. The initial region has a slightly anisotropic shape, 
but this is quickly counteracted by the shear (cf. Figs. 2 and 3). The density contours are 
pulled into a strongly triaxial shape, as one axis is pushed in by the shear and a second axis 
is pulled out; the third axis is somewhere in between. As the central peak collapses along 
its short axis, which is more likely to be considered as part of the object: the high density 
regions located far out on the long axis or the lower density regions located nearby on the 
short axis? For the formation of a galactic halo, the nearby regions are more relevant, since 
they are positioned closer to the center of mass and therefore undergo shell crossing earlier. 
Moreover, the sections of the high initial density contours out on the long axis are susceptible 
to being separated from the halo by fragmentation or other instabilities (Merritt & Hernquist 
1991). Thus, the initial isodensity contours around the peak may not be appropriate tracers 
of the volume that eventually makes up the halo. 

An additional problem is that while a galaxy corresponds to the region surrounding a 
high-density peak, the peak itself collapses at high redshifts and the surrounding matter 
accretes onto this seed. Hence, the properties of the collapsed peak do not reflect the 
properties of the final halo. In a spherical collapse model, one hopes to remedy this problem 
by picking the radius of the object so that the overdensity inside that radius corresponds to 
a collapse at present. However, when the effects of shear and non-sphericity are included, 
this procedure is no longer valid. 

For the above reasons, we feel that the ellipsoid model does not provide a fully sat- 
isfactory description of the dynamics of typical galactic halos. However, the model does 
substantially better in describing low-spin objects at high redshifts (cf. Paper II). Such ob- 
jects are generally located in regions with low shear, so that their collapse is more spherical, 
as demonstrated in Figure 8. This fact tends to reduce the concerns about the proper initial 
volume for the system. The population of low-spin objects is of particular interest as po- 
tential environments that favor the formation of massive black holes, since in these systems 
the hydrodynamic collapse of the gas is not inhibited by the centrifugal barrier as in typical 
systems (e.g., typical disk galaxies are larger than their Schwarzschild radius by 6-8 orders of 
magnitude because of rotational support). In Paper II we show that the existence of low-spin 
systems at high redshifts can in principle account for the seeds of quasar black holes with 
masses ^ 10 6 Mq and a comoving density of bright galaxies. Appendix A then shows that 
if a black hole forms in the initial collapse of a high-a peak on the ~ 1O 8 M mass scale, it 
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is likely to be surrounded by a highly overdense region even on the mass scale of a galactic 
bulge. Due to the proximity of the centers of mass of the two systems, the black hole will 
sink by dynamical friction to the center of the bulge system after it forms. The later collapse 
of the surrounding region would fuel the black hole and result in the quasar activity (Loeb 
& Rasio 1994). In the application of the ellipsoid model to the formation of quasar black 
holes, we are free to consider only the inner part of the collapse and to neglect subsequent 
accretion, since the formation of the seed black hole occurs prior to this accretion. The main 
remaining uncertainty in this treatment is the omission of higher multipole torque couplings. 

Given a primordial power-spectrum of Gaussian density perturbations, the ellipsoid 
model makes definite predictions about the statistics of shapes and angular momenta of 
overdense regions in the universe. It would first be useful to compare the quantitative pre- 
dictions of our model to the statistical properties of collapsing regions in high-resolution 
simulations (Bertschinger 1993 and references therein). Such a comparison would require 
large simulated volumes (<; [lOOMpc] 3 ) in order to obtain reasonable statistical samples of 
rare objects that collapse today (M ^ 1O 15 M ). The prominence of filaments and sheets 
relative to quasi-spherical structures could also be tested observationally using galaxy sur- 
veys (e.g., Geller & Huchra 1989; Maddox et al. 1990; Saunders et al. 1991; Shectman et 
al. 1992; Strauss et al. 1992; or the future DSS survey, Gunn & Knapp 1993) to infer the 
smoothed mass distribution on scales larger than the virialized cores of clusters of galaxies. 
A non-spherical analysis of this type can provide constraints which are complementary to 
the conventional results obtained by applying spherical filters to galaxy surveys. 



We thank John Dubinski and David Weinberg for useful discussions. D.J.E. was supported 
in part by a National Science Foundation Graduate Research Fellowship. 
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APPENDIX A: PROBABILITY DISTRIBUTION FOR THE INITIAL DATA 

Our model has 6N + 5 real numbers as its initial data, each of which is defined as some 
integral over the initial density field. In this Appendix we derive the probability distribution 
for these initial data in terms of the properties of the Gaussian random field of initial density 
perturbations. 

We denote the radii of the boundaries between the shells as Rq, R\, . . . , Rn, where Rn = 
oo. The values to be derived from the field are (m = 0, ±1, ±2): 

* {Rn) = {±^Rj) I d3r5{r) forn = °. 1 .---.^- 1 ; ( A1 ) 

\r\<R n 

q2m = Pb J d 3 r5(r)r 2 Y£ m ; (A2) 
|r|<.Ro 

42 = P6 J d 3 sY 2 * m 5(s)s- 3 forn = 0,l,...,iV-l. (A3) 

7?„<|s|<oo 

/From these values, we can easily transform back to our model input variables in equation 
(34) and (33), noting 

4S = 4r 1} -4S forn = l,2,. (A4) 

Let us first switch to a multipole expansion of the Gaussian random field, as described by 
Binney & Quinn (1991). We replace the random field 8{x) by a new set of random functions 

{<W*;)}: 

nr oo i °9 

5 W = V"EE / dk8 lm {k)ku{kr)YU0A), (A5) 

where jl are the spherical Bessel functions and (r, 9, </>) are spherical coordinates. Binney & 
Quinn (1991) show that the set of functions {<W} is Gaussian distributed as 



2P(k) 



(A6) 



where P(k) is the power spectrum. This means that the functions 5g m (k) are independent 
and that all of them have the same simple probability distribution. 
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When we insert equation (A5) into equations (A1)-(A3), the angular integrals may be 
easily done. The q2m and coefficients depend only on (W(^), and 5(R n ) depends only 
on doo(k). All the other £ in the expansion do not appear, so we do not need to determine 
those 8t m {k). Furthermore, the values we seek fall into six independent sets, since none of 
the integrals in equations (Al)-(A3) mix different {l,m) sets. 

We next perform the radial integrals by using the identities (Abramowitz & Stegun 1957) 



We thus find 



Similarly, 



j- z [z n+1 Uz)] =z n+1 j n - 1 {z), 



, — oo R n 

S(R n ) = y ^ J dk J drSo Q (k) r 2 kj (kr) 

o o 

oo kRo 



AixRl 



\J^- J dkk 2 5 o(k) J dzz 2 j (z) 





oo 



dkk^Sooik^kRnfhikRn) 



(AT) 



(A8) 



47Ti?3 



oo 

3 f „ x , Ji(kR n ) 



i— 00 

Q2 m = Ph Rl\Jl J dk5 2m (k) 3 -^^- (A9) 



and 

— 00 

b^ = Pb f-Jdk5 2m {k) J -^^. (AlO) 
o 

Each desired quantity has been written as an integral over the Si m (k), a property that will 
be used later in finding the relevant probability distribution. 

Before developing the joint probability distribution for the full problem, we wish to illus- 
trate the method by working a simpler example first. We consider the average overdensities 
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inside two different radii and ask: given a value for the overdensity inside the smaller radius, 
what is the distribution for that at the larger radius? Denoting the radii R\ and R2, we 
have from equation (A9): 

00 

7i = 4J*M*) il(Jkfll) 



( A11 ) 

3 [..^uJlW 



82 = — ^= y dk5oo(k)- 



R2 



Because the functions 8t m {k) are independent and because we need only consider one at 
a time, we will suppress the £m label. The form of the integral in equation (A6) prompts us 
to consider 6(k) as a vector in a function space and to define the inner product as 



00 





for two arbitrary functions a and b. Next, define: 



/dk 
— a*(k)b(k), (A12) 



(AlH) 



Then we have 

3 



, - -(Qil^oo) (A14) 
ny 1 

and similarly for $2; note that the P(k) in the definition of Q cancels the P(k) in the inner 
product. 

Now let us consider some basis n,j(k) that is orthonormal under our inner product. Then 
we may write S(k) as a linear combination of these basis vectors, 

5(A0 = I>™#)- (A15) 
3 

Substituting this form into equation (A6), we find the probability for a given set of f3j to 
occur, 

P^ocexp^-i^^^ln^j =exp^-i^|/? J fj = ]J exp (-^|^| 2 ) ■ (A16) 

This means that each f3 is independently distributed with a uniformly distributed phase and 
a Gaussian distributed magnitude. For an m = function as we have here, (5 is real and 



29 



is simply a Gaussian deviate, namely it is drawn from a Gaussian distribution with a zero 
mean and a unit variance. 

We would like to express Q\ and Q2 in terms of this ortho normal basis. Since we have 
not yet constructed the basis, we choose it to be simple. The first vector n\ will be Qi 
normalized, 

ni = -. Ql . (A17) 
Then, we pick 112 to be the normalized component of Q2 orthogonal to ni, which is 

n 2 = , g2 ' (ni|g2)ni (A18) 



\/ (Q2\Q2) - (ni\Q 2 )' 



We may then complete the basis any way we wish; none of the other vectors will include any 
component of Q\ or Q2. Inverting (A17) and (A18), we find 



Qi = V {Qi\Qi)ni, 

I 2 (A19) 

Q 2 = (ni\Q2) n\ + y (Q2IQ2) - {n\\Q2) n 2 . 

Putting equations (A9) and (A15) into (A14) and using the orthonormality of the basis, 
we find 

3 



Si = —j=y/{Qi\Qi)Pi (A20) 

7TV 2 



and 



(A21) 



62 = ^7! ( (ni|g2) Pl + V( < 52|Q2)-<ni|Q2) 2 /9 2 
= ^V(Q2W) (tA + y/i - 1^2) , 

where 7 is the dimensionless overlap defined as 

To determine 5\ and ^2 randomly, we simply find two Gaussian deviates (3\ and P2 and 
combine them as indicated. Alternatively, if we are given S\ and asked to determine 82 given 
this constraint, we fix the value of f3\ to produce 8\. Despite this constraint, ^2 remains a 
Gaussian deviate since it is independent of f3\, and so we may take a random value of f3 2 
along with our fixed value of (5\ to find the constrained distribution of 8 2 - 



30 



Equation (A20) can be expanded to give 



<JXJ o 



{6l) = 2^J dkP(k) llh ' (A23) 



o 



the expectation value of which is the standard expression for {(5M/M) 2 ) for a spherical 
top-hat window function of radius R\ [cf. Eq. (36)]. Denoting this expectation value as 
(oi) 2 and similarly for (02) 2 , we find v \ = = A and z/2 = ^2/o"2 = 7/3i + a/1 — 7 2 /?2, 

so a high-cr peak on one mass scale will correspond to a high-cr peak on the other mass scale 
if 7 is close to I. 

For the f2 = 1 CDM model and mass scales of 1O 8 M and 1O 1O M , we calculate 
7 Ki 0.74. This means, for example, that given a 3cr peak at 1O 8 M , the surrounding 
10 10 M region has an overdensity of z/2 = 2.22 ± 0.67(3, where (3 is a Gaussian deviate. This 
result has interesting physical implications to the problem of the origin of quasar black holes 
discussed in Paper II. If a black hole forms in the initial collapse of a high-cr peak on the 
~ 10 8 M mass scale, it is likely to be surrounded by a highly overdense region even on the 
mass scale of a galactic bulge. 

We now need to generalize the above method to an arbitrary number of vectors. We 
seek p values from the field, denoted fi, f2, ■ ■ ■ , fp, and each may be written as some integral 
over a particular 5i m (k); only one (£, m) is involved and so we suppress the index here as 
well. Then, we define the Qj{k) by the relations 

00 

fi = j ^S(k)Qj(k), (A24) 


where all of the Qj are real-valued functions; and define the matrix M, 

00 

/dk 
p^Qi(k)Qj(k). (A25) 



We now seek an ortho normal basis {nj} with the property that 



p 



Qi = Y, L i3 n h (A26) 
where L is a lower triangular matrix and the remaining nj are not used. Substituting 
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equation (A26) into (A25) and using the ort honor mality property, we find 




(A27) 



Hence, our basis must be defined by M = LL 1 . Can we find such an L? In fact we can, for 
if M is postive-definite symmetric matrix, this is the Cholesky decomposition (Press et al. 
1992, §2.9). M is obviously symmetric and must be positive-definite if P(k) > for all k, 
and so we may find a unique L. 

Next, we take {nj} as the basis for S(k). Then, as before, we decompose by equation 
(A15) and reach the probability distribution in (A16). Substituting this result, as well as 
equation (A26), into (A24), we get 



This completes the solution. One needs only to specify p Gaussian deviates and combine 
them in the indicated way in order to get realizations of the p values fj. To construct the L 
matrix, one must perform p(p + l)/2 one-dimensional integrals as given in equation (A25). 
This task can be made even simpler by doing the integrals together since they have many 
functional evaluations in common. 

In the particular case of the initial conditions for the ellipsoid model, we have six sets 
of parameters to determine, one for £ = and five for £ = 2. But the five £ = 2 sets are 
equivalent and will share one L matrix. Also, the integral forms for £ = in equation (A8) 
are identical to those for £ = 2 as given in equation (All). The only difference between 
the two sets is the one additional initial condition presented in equation (A9). In order 
to compute L for the £ = set, we need not calculate any new integrals; we simply take 
the required integrals as a submatrix of M from the £ = 2 calculation, and compute a new 
Cholesky decomposition. 

Because the original field <5(x) was real, the complex functions 5g m (k) have the usual 
spherical harmonic relation between them, i.e. 5g m (k) = (— l) m 5p _ m (k). Thus, for m ^ 
there are four real functions: the real and imaginary parts of 5i m {k) and 5£- m (k), only two 
of which are independent. Returning to equations (14) and (20), where £ = 2, we choose to 

work with the real and imaginary parts of m — 1 and m = 2 as well as the real term m = 

(n) 

as the five independent real parameters. This means that instead of asking for q2 m and 



•DO 




(A28) 
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as in equations (A2) and (A3), we compute Reg22, bmg22, Re 921, Img2i, and q2o, and the 
same for ti- n \ This is easy to do by simply taking the quantities Re 622(h), Im^22, as 
the independent random functions. The t = 2 values again fall into five independent sets, 
each being treated as above. There is however one modification, as seen from the original 
multipole field probability distribution in equation (A6): 



P[5(x)] oc exp 



00 

£,m q 



\hm\ 2 

2P(k) 



exp 



2(Refe) 2 + 2(Im 



2(Re^_ 1) ) 2 + ... + (^ ) 



2P(k) 



(A29) 



This indicates that the real and imaginary parts of the m/0 terms actually are distributed 
with half the variance of the m = term. The simplest way to incorporate this fact is to 
compute all the £ = 2 sets as for the m = case and then take the [3j to have variance of 
one-half rather than unity for m^fl. 

Finally, there is a choice as to what order to place our fj in the vector. Since the matrix 
L is a lower triangular, the value fi requires only one Gaussian deviate to be determined. 
As discussed in §4, there are additional constraints that we place on the initial conditions. 
Since these are imposed simply by rejecting any set of initial data that do not satisfy them, 
it is most efficient to pick the most discriminating test to be f\. In particular, the test that 
<5(i?o) > v m i n <j is very stringent and we arrange the t = functions so that this condition is 
the first basis function. This way, we need only find one Gaussian deviate before applying 
this test. 
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APPENDIX B: TREATMENT OF AN AXIS COLLAPSE 



The inclusion of a collapsed axis is complicated by the facts that the ellipsoid is rotating 
and that the columns of the matrix A need not be orthogonal. When we halt the collapse 
of an axis, we wish only to end the radial collapse and to leave the tangential velocity 
unchanged. 

First, we must identify that an axis has reached 40% of its turn-around value. At each 
time step, we diagonalize AA T = QAQ T to find the axis lengths, which are in y/X. We put 
these in ascending order and compare each to the previous maximum length for the short, 
middle, and long axes. If the new length is longer, we save it as the new maximum. If 
the new length is less than 40% of the maximum for that axis, we assume that the axis has 
collapsed and save the column number so as to be able to refer to it in Q and A. This method 
of identifying the axes by sorting the lengths works well because the axes are well-separated 
soon after the beginning of the integration and well before turn-around. 

If a particular axis has collapsed, we remove the radial component of the velocity in 
that direction and we alter the right hand side of the equation of motion (8) so that there 
is no inward force in that direction. The steps for each alteration are equivalent, so we only 
describe the first in detail. 

Suppose we begin from the diagonalization AA T = QAQ T and want to fix the a-th axis. 
The position of mass elements within the ellipsoid may be written as r = Q\/Xs, where s 
is a vector on the unit sphere. Then the frozen axis is s || e Q , the Cartesian basis vector. 
The velocity field is v = AA~ 1 Q\/Xs, where A = dA/dt. We now rotate this velocity to the 
principal axis frame of the ellipsoid, which requires left-multiplying by Q T . This is obtained 
by defining 



where Vjk is the jth component of the velocity at the point in the direction of the kth axis, 
and the axes are numbered according to the columns of Q. 

Now we can remove the a-th component of velocity along the a-th axis by removing the 
aa component of V. We therefore define a new matrix W such that 



v = q t aa~ 1 qVX, 



(Bl) 




if (3 — a and 7 = a, 
Vg 7 otherwise. 



(B2) 



Then we need to connect this matrix back to the new velocity matrix A new . We do this by 
solving equation (Bl) for A after replacing V by U . Then we replace the old matrix for A 
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with 



A new = QUVX Q T A. (B3) 

If more than one axis has collapsed, we modify the definition of U to eliminate all diagonal 
elements corresponding to the collapsed axes. This algorithm works equally well for the forces 
by substituting A = d 2 A/dt 2 for A everywhere. Since this procedure does not modify the 
tangential velocities, it preserves angular momentum; however it does change the vorticity. 



APPENDIX C: RELATION BETWEEN SHEAR AND OVERDENSITY 

In this Appendix we present a calculation separate from the dynamical ellipsoid model, 
but yet based upon the formalism discussed in this paper. As before, we begin by dividing 
the universe into two pieces through a spherical boundary of radius R around the origin. 
We now wish to compare the statistical properties of the average overdensity in the interior 
region 5(R) [cf. Eq. (15)] to those of the quadrupole shear resulting from the exterior region 
a2m [cf- Eq- (2)]- We will show that these quantities are drawn from independent Gaussian 
distributions whose variances are related by a constant of proportionality that is independent 
of the power spectrum or mass scale. 

iFrom equations (Al), (A3), (A8), and (All), we can easily relate 5 and a2m to integrals 
over the coefficients of the multipole expansion of the Gaussian random field: 

oo 

5 = ^Jdk5 00 (k) J ^ (CI) 
o 

and 

I — 00 

a 2m = p b ^l J dk5 2 m(k) 3 -^^. (C2) 
o 

Because each of these six quantities depends upon different Si m , they will be independent. 
By the methods of Appendix A, we find that the six quantities are Gaussian distributed 



35 



with a zero mean and the variances 



6) = { W=2) J dkP ^{^)> (C3) 



2 

<(a 20 ) 2 > = Lyf) jdkpw^hmy, 



(C4) 



and 



<(Rea 22 ) 2 > = <(Ima 22 ) 2 > = <(Rea 21 ) 2 ) = <(Ima 21 ) 2 ) = 1 <(a 20 ) 2 > . (C5) 



The last equation follows from equation (A29) and the discussion around it. 

The key result of this derivation is that the integrals in equations (C3) and (C4) are 
identical. We can therefore pull all of the dependence on the power spectrum into one 
constant. We define a as the r.m.s. amplitude of 5M/M according to equation (36). We may 
then write 5 = ais, a 2 o = (2y/jfpi,/3)<jzi, Rea 22 = (v / 27rp&/3) crz 2 , Ima 22 = (\ / 2%pf ) /3)az3, 
Rea 2 i = (\^2npi ) /3)az4, and Ima 2 i = (v / 27rp6/3)<T-25, where v and the {zj} are all Gaussian 
deviates of unit variance. 

We next consider how the inner region would evolve given these values as the initial 
conditions. In particular, we apply the Zel'dovich approximation to find whether all the 
material of the interior region will turn-around and form a spatially bounded object or 
whether the region will expand forever in at least one direction. For an Q = 1 universe 
and under the Zel'dovich approximation, this question is simply a matter of whether the 
initial peculiar velocities have inward radial components or not. This, in turn, depends on 
the initial gravitational potential (as we assume growing mode velocities), which we may 
construct from equations (12) and (14) to get the matrix 

^ = ^^^IaP + Sap), (C6) 

where / is the identity matrix and 

5 | ^ -z 2 -f 3 z 5 . (C7) 

We want $ to have positive eigenvalues, which means that the eigenvalues of S must be 
greater than —v. The probability that a spherical region of overdensity v will expand without 
bound in one direction due to shear depends only on the distribution of the random matrix 
S. All reference to the power spectrum has been removed. 
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Table 1 shows the cumulative distribution function for the most negative eigenvalue 
—[/, (/j, > 0) of the matrix S. This quantity acts as a negative u, counteracting the actual 
v = S /a that one normally identifies as the parameter controlling gravitational collapse. The 
importance of shear through the collapse is evident from the prevalence of large values of \i. 
For low Q universes, one might alter the peculiar velocity requirement as discussed in §4 [cf. 
Eq. (35)]; this leads to the requirement that v — (j, > 5 cr it/cr. 

The above calculation provides insight to the acceptance rate for the velocity constraint 
that we use in the ellipsoid model. In this model, there is an additional term in the potential 
coming from the anisotropy of the ellipsoid. This term tends to slightly increase the accep- 
tance (i.e. resist the shear) because the long axis of the ellipsoid tends to be aligned with the 
most negative shear axis. The above calculation also has implications to the Press-Schechter 
formalism (Press & Schechter 1974), in which spherical regions with a sufficiently high v are 
all assumed to turn into objects at the appropriate collapse time. However, as it stands the 
model can only eliminate objects on the basis of shear, some of which would presumably be 
re-introduced through mergers of lower mass systems or fragmentation of unbound filaments. 
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FIGURE CAPTIONS 

Fig. 1: Probability distribution of the spin parameter P(X) for the mass scales of 10 15 M 
and 10 12 M and u m i n = 2.5 in an f2 = 1 universe. Each histogram includes 2 x 10 5 random 
realizations of the initial conditions. The average values (A) = J Q °° XP(X)d\ are also shown. 

Fig. 2: Distribution of axis ratios for an ensemble of 5000 systems of mass 1O 15 M and 
v = 2 in an Q = 1 universe. Panel (a) shows the initial distribution and panels (b)-(d) 
show the resulting evolution at particular times. Panel (e) shows the distribution at time 
^max, the end time of the integration defined by the collapse time of a homogeneous spherical 
perturbation with the same initial density as the ellipsoids. Panel (f) shows the mean values 
of the distributions in the other panels. We define the axis lengths as a > b > c. To find the 
redshifts we assumed a bias parameter of 1.3 and Hq = 50km/s/Mpc. 

Fig. 3: As in Figure 2, but with the initial ellipsoid replaced by a sphere of the same density. 
The external shear was unchanged. 

Fig. 4: As in Figure 2, but for 1O 12 M regions. 

Fig. 5: As in Figure 2, but for an open universe with = 0.2, Ho = 50km/s/Mpc, and no 
bias (b = 1.0). For 1O 15 M ellipsoids, we pick v = 2.2 in order to match the final collapse 
redshift of Figure 2 and then pick the time slices to match the other redshifts of that figure. 

Fig. 6: As in Figure 5, but after removing the constraint on the initial peculiar velocities. 
This allows the inclusion of objects with at least one axis that is likely to expand without 
bound. 

Fig. 7: As in Figure 2, but for 1O 8 M regions. 

Fig. 8: As in Figure 7, but showing only the 232 low-spin objects from an overall sample 
of 2 x 10 5 systems. 

Fig. 9: Four individual examples of ellipsoids evolved by the model [panel sets (a)-(d)]. The 
left panel of each set shows the lengths of the three semi-axes as functions of time. The middle 
panel of each set shows the overdensity of the ellipsoid (solid line) as a function of redshift 
z. Also shown are the overdensity of a spherical top-hat perturbation of equal initial density 
(dotted line) and the overdensity predicted from applying the Zel'dovich approximation to 
the initial potential (dashed line). The right panel shows the spin parameter A as a function 
of time. 
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/i P(fJ, < /io) 

.20 .01% 
.31 .1% 
.53 1% 
.91 10% 
1.19 25% 



/i < no) 

1.55 50% 

1.96 75% 

2.36 90% 

3.11 99% 

3.69 99.9% 



Table 1: The probability that the most negative eigenvalue —\x of the matrix S in Appendix 
C is greater than a particular threshold — /io- 
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